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Abstract 

We have carried out simulations of the coalescence between two relativistic clus- 
ters of collisionless particles using a 3D numerical relativity code. We have 
adopted a new spatial gauge condition obtained by slightly modifying the mini- 
mum distortion gauge condition proposed by Smarr and York and resulting in a 
simpler equation for the shift vector. Using this gauge condition, we have per- 
formed several simulations of the merger between two identical clusters in which 
we have varied the compaction, the type of internal motion in the clusters, and 
the magnitude of the orbital velocity. As a result of the coalescence, either a new 
rotating cluster or a black hole is formed. In the case in which a black hole is not 
formed, simulations could be carried out for a time much longer than the dynam- 
ical time scale, and the resulting gravitational waveforms were calculated fairly 
accurately: In these cases, the amplitude of gravitational waves emitted can be 
~ 10~ 18 (M/10 6 M Q ) at a distance 4000Mpc, and ~ 0.5% of the rest mass energy 
may be dissipated by the gravitational wave emission in the final phase of the 
merger. These results confirm that the new spatial gauge condition is promising 
in many problems at least up to the formation of black holes. In the case in which 
a black hole is formed, on the other hand, the gauge condition seems to be less 
adequate, but we suggest a strategy to improve it in this case. All of the results 
obtained confirm the robustness of our formulation and the ability of our code 
for stable evolution of strong gravitational fields of compact binaries. 



§1. Introduction 

The coalescences of relativistic binaries, are the most promising sources f on kilometer- 
size laser interferometers such as LIGO^'i VIRGO, & GEO,& and TAMA,!) which 
will be in operation in the next five years, as well as for future space-based interfer- 
ometric gravitational wave detectors such as LISA.B 1 When a signal of gravitational 
waves from such compact objects is detected, it will be analyzed using matched filter 
techniques, and a variety of astrophysical information will be extracted.^ In order 
to apply this technique, theoretical templates of gravitational waveforms will be 
needed, and this motivates the recent theoretical study of gravitational waves emit- 
ted from coalescing binaries. The post-Newtonian approximation has been shown to 
be a powerful tool for the study of the inspiraling phase, and a large effort devoted 
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to this study has produced many satisfactory results. IIP However, in order to study 
the final phase of the coalescence, in which strong general relativistic effects domi- 
nate, no approximation of general relativity is valid, and 3D numerical relativistic 
simulations appear to be the only promising approach. 

Although a few results have been obtained recently, B> there are still a number 
of fundamental questions that need to be addressed. Among them, it is particularly 
important and urgent to find appropriate gauge conditions which allow one to per- 
form stable simulations and to extract gravitational waves accurately. As shown by 
many numerical simulations, the maximal slice condition seems to be an adequate 
time slice condition at least up to the formation of black holes. However, a similarly 
reliable spatial gauge condition in 3D numerical relativity has not yet been studied 
sufficiently More than 20 years ago, Smarr and YorkEP proposed a minimum distor- 
tion (MD) gauge condition, which is expected to have several excellent mathematical 
properties. From a numerical point of view, however, the MD gauge condition re- 
quires the solution of a complicated vector elliptic equation, and this can be a very 
time-consuming task. 

Their MD gauge condition was contrived from the idea that the global change 
rate of the conformal spatial metric should be minimized.™ This is a rather restrictive 
condition, and we believe that a gauge condition in which the global change rate 
is not exactly minimized but instead only approximately minimized can prove to 
be appropriate, as long as the geometrical distortion of the three space does not 
increase quickly. Thus, in this paper, we propose an approximate minimum distortion 
(AMD) gauge condition which is similar to Smarr- York's MD gauge condition and 
in which the global change rate is expected to be only approximately minimized. In 
a numerical simulation, the AMD gauge condition is much more easily imposed than 
the original MD gauge condition. 

In a previous paper, 113 we presented results from numerical simulations of black 
hole formation in which collisionless particles were used as matter sources of the 
Einstein equation. That study demonstrated that our formulation is robust, that 
it allows for stable numerical evolution, and that it can be applied to problems 
such as the formation of black holes from collision of two clusters, or gravitational 
collapse. In this paper, we use our code to simulate the merger between two identical 
binary clusters, adopting an AMD gauge condition, and we demonstrate that this 
formulation with a new spatial gauge condition allows for stable numerical evolutions 
of merging binaries and for fairly accurate calculation of gravitational waves. 

Numerical simulations of the merger between relativistic clusters may become a 
source of important information in astrophysics and gravitational wave astronomy. 
It has been in fact proposed that supermassive black holes in galactic centers might 
be formed through complex phenomena such as gravitational collapse and bar insta- 
bility, involving the coalescence between clusters of collisionless compact stars. liiP If 
such phenomena occurred frequently in the early universe, they could have produced 
low frequency gravitational waves which would be detected by the planned interfer- 
ometric detectors in space. EP In particular, the coalescence of two compact clusters 
could be an important source. The simulation presented in this paper may play a 
role in predicting the amplitude and frequency of the gravitational waves emitted 
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during such coalescences. 

The paper is organized as follows. In §2, we present the basic equations to 
be solved in the numerical simulation. In §3, we describe the spatial AMD gauge 
condition adopted in this paper and briefly discuss its properties. In §4, we describe 
the methods employed in order to analyze the gravitational waves radiated. In §5, 
after we briefly describe the method to set the initial conditions, we present numerical 
results for simulations of the coalescence of two clusters. Section 6 is devoted to 
summary. Throughout this paper, we adopt units in which G = 1 = c. Latin and 
Greek indices denote spatial components (1 — 3) and spacetime components (0 — 3), 
respectively. As spatial coordinates, we use the Cartesian coordinates x k = (x, y, z) 
with r = \J x 2 + y 2 + z 1 . 



§2. Basic equations 



JDur formulation of the Einstein equation is described in detail in previous papers. 
E3)'I12P We here briefly recall the basic equations and refer the reader to the above 
cited works for details. 

We write the line element in the form 

ds 2 = (-a 2 + (3 k l3 k )dt 2 + 2(i i dx i dt + lij dx i dx j , (2-1) 

where a, (3 l (Pi = 'jijfl 3 ) and 7^ are the lapse function, shift vector, and 3D spatial 
metric, respectively. We define the quantities 

7 = det( 7 tf) = e 12 t (2-2) 
jij = e'^jij, i.e., det(Tjj) = 1, (2-3) 

Aj = e-^(K i:j - \lijK\ (2-4) 

where Kij is the extrinsic curvature, and K = K k k . The indices of Aij and/or are 
raised and lowered in terms of 7^ and 7*- ? . Numerical computations are carried out 
using 7jj, Aij, (ft and K as dynamical variables rather than 7^' and Kij. Hereafter, we 
use Di and Di as the covariant derivatives with respect to jij and 7^, respectively. 

Since our matter source is represented by collisionless particles, its energy mo- 
mentum tensor is written as 




where m p is the rest mass of each particle (we assume that all particles have the 
same rest mass), x J a denotes the position of the a-th particle, N is the total particle 
number, is the four-velocity of a particle, and 5^(x^ — x J a ) is the Dirac delta 
function in the 3D spatial hypersurface. Note that the rest mass density and the 
total rest mass (i.e., conserved mass) M* are written as 

N 

P* = ™ P 5> 3 V-4), (2-6) 

0=1 

M* = 4A^m p , (2-7) 
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where the factor 4 in Eq. (2-7) arises because in this paper we assume a plane sym- 
metry with respect to the equatorial plane (z = 0) as well as a 7r-rotation symmetry 
around the z-axis, and treat only a quadrant region (see below). 

The equations of motion for the particles are derived from the geodesic equations 
and are written in the form 



dui 



UjU k j k 

~2^ 7 >* 



where Qj = djQ = dQ/dxK Once itj is obtained, u° is determined from the nor- 
malization relation of the four-velocity as 



(era ) 2 = 1 + j t3 UiUj 



(2-9) 



We also have a relation between a coordinate velocity dx l /dt = vf/iP and Uj as 



dx l 
~dt 



(2-10) 



Equations (2-8) and ( 2TC| ) are the basic equations for evolution of the particle's 
position and velocity. 

As customary in a 3+1 decomposition, the Einstein equation is split into the 
constraint and evolution equations. The Hamiltonian and momentum constraint 
equations are given by 



R - AijA'i + -K 2 = IQttE, 



DiA 1 



-DjK = 8vrJ,-, 



where 



N 



E = m p Y J {u°)aae- & n^(x k -x k a ), 

0=1 

N 

J i = m p J2(ui)ae- 6<l, 5 i3) (x k -x k a ), 



(2-11) 
(2-12) 

(2-13) 
(2-14) 



a=l 



and R is the scalar curvature with respect to jij. We solve these constraint equations 
only when setting initial conditions. 

The evolution equations for the geometric variables 7^, eft, A^j and K are given 
by0> 



(d t - I3 k d k )% = -2aA i:j + + % k p k ti - -7ii/3 fc fe , 



(2-15) 



(d t - k d k )A ij 



-4<j> 



a Ri 



+a(KA lJ - 2A ik A k ) + ^A kj + p k tj A ki - ~p k k A 



ij 
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-Snae-^^-^S^, (2-16) 

(ft - f3 k d k )(P = 1 (-aK + f3 k ^ , (2-17) 
ikn_\js ( A. Aii i i n. r>fc„, i ^„,/p i c ^ 



(ft - /? fc ft)ir = aMyA* + -2^ J - D k D k a + 47ra(£ + (2-18) 
where Rij is the Ricci tensor with respect to "fij, and 



We also solve the following equation for Fj = 5 J 7« 



(ft - /3 fc ft)F i= 2a(/^i ifej - + / fc ^.I ife - ^i% M + 60, fc A fc i - ^K, 

-2<5**a, A + <F^% fc + (7^ + 7^, 4 - |%-/9 , ,i) I fc^'* 
-16vraJi . (2-20) 

Here <5jj denotes the Krpnecker delta, h{j = jij — Sij, and fy = ^f 3 — 5 iJ . Fi is used 
to compute Rij and R.i3 As argued in previous papers, 113' c3 the introduction of 
Fi is necessary for stable numerical computation. 

As the time slice condition, we use an approximate maximal slice condition. 
Namely, we determine a by demanding that the right-hand side of Eq. (2- 18L be- 
comes approximately zero, using a strategy discussed in a previous paper. i3 A 
discussion of the spatial gauge condition used in this paper will be presented in the 
next section. 

The numerical methods employed in solving Eqs. ( p -8| ) , ( 2TC| ), (2-15)-(2T8) and 



(2-20), and for finding, the apparent horizon are almost the same as those discussed 
in a previous paper, 113 except for the following two differences. The first difference 
is due to the appearance of the transport terms such as —fftdklij hi the evolution 
equations of the geometric variables. This is because we choose a non-zero (5 k (in 
a previous paper EU 1 we adopted f3 k = 0). A discussion of our strategy for handling 
such transport terms is described in Appendix A. The second difference is due to the 
different underlying symmetries assumed here. In this paper, we solve equations in 
the quadrant region —L < x < L, < y, z < L, where L denotes the location of the 
outer boundaries, assuming 7r-rotation symmetry around the z-axis as well as the 
plane symmetry with respect to the z = plane (in a previous paperlilP we assumed 
triplane symmetry). We note that we perform simulations of binary clusters whose 
mass centers orbit in the z = plane (see §5). Boundary conditions in the y = 
plane are as follows: 

Q(x,0,z) = Q{-x,0,z), (2-21) 
Q A (x, 0, z) = -Q A (-x, 0, z), Q A (x, 0, z) = -Q A (-x, 0, z), (2-22) 
Q z (x,0,z) =Q z (-x,0,z), Q z (x,0,z) =Q z (-x,0,z), (2-23) 
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Qab{x, 0, z) = Qab(~x, 0, z), (2-24) 
Q Az (x, 0, z) = -Q A z{-x, 0, z), (2-25) 
Q zz (x,0,z) =Q zz (-x,0,z), (2-26) 

where A = x or y, and Q, Q J (or Qj) and Qij denote arbitrary scalar, vector and tensor 
quantities, respectively. Note that the boundary conditions at the outer boundaries 
are the same as those in a previous paper t3 except for that of Fi for which we 
impose Fi = 0(r~ 3 ) in this paper. It is also convenient to define the quantities 

N 

j(r) = m p ^2(xu y - yu x ) a for r a < r, (2-27) 

a=l 

m*(r) =m p N p (r), (2-28) 

where r a = x\ + y'^ + z* and N p (r) denotes the number of particles which are in 
a radius r. We consider j(r) to be the approximate z-component of the angular 
momentum within r, while m*(r) denotes the total rest mass within r. 

§3. Spatial gauge condition 

The MD gauge condition proposed by Smarr and Yorki) can be written as 

D\e^da i3 ) = 0, (3-1) 

or, equivalently, as 

D^e^dt^j) = 0. (3-2) 
More explicitly, their MD gauge condition reduces to an equation for (3 k : 

... 4 . 

-2AijD % a- -aDjK = WnaJj. (3-3) 
3 

Here = "fki/3 l and f3 k = f3 k . Smarr- York's MD gauge condition is also derived by 
minimizing the following action I with respect to [3 k on three spacelike hyper surf aces: 



1= / d 3 x(da i3 )(da kl )f k f l e 6 t (3-4) 



Thus, by choosing the Smarr- York's MD gauge condition, the global change rate of 
jij based on the action I is minimized in every three hypersurface. 

Since I does not have a special physical meaning, we may consider an alternative 
gauge condition by slightly changing the definition of the action. In particular, we 
define another action V as 

I'= [ a*x{da l3 ){da kl )f k y l . (3-5) 
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This corresponds to defining the action in the conformal three space in which the 
determinant of the metric det(7y) is unity. In this case, by taking the variation of 
I 1 with respect to /3 k , we obtain a different MD gauge condition, 



0, 



(3-6) 



or more explicitly, 



D*[D i (3 j + D j /3 i --%D k p 



-aDjK = WnaJj. 
3 



(3-7) 



A merit of this condition may be that we do not have a coupling term between (3 k 
and 4>, and hence the equation for (3 k is slightly simplified compared with Eq. (3-3). 
However, it is still complicated to solve in numerical computation. 

Although it is desirable to adopt a shift vector in which / or I' is exactly min- 
imized, we believe that we do not always have to use such a shift vector. We may 
probably use another shift vector in which / or /' is approximately minimized. Using 
this idea, we derive a gauge condition by slightly changing Eq. (3-7): We rewrite the 
covariant derivative operator acting on (3 k as a partial derivative; i.e., we solve 
the following equation to determine /3 fc (= f3 k ): 



5. 



l3 Ap + 



2A ij (D i a-6aD i 



-aDjK = 167ra J 



(3-8) 



Here A is the Laplacian in the flat space. In this case, Eq. ( |3-8[ ) can be rewritten 
into simple elliptic equations for a vector P, and a scalar r/ using the transformation 



where Pj and r] satisfy 



Arj = 



S. 



it 
Six 



and 



Si = 16vrQJi + 2A i0 X&a - 6aD j 



4 ~ 
+ -aDiK. 
3 



(3-9) 

(3-10) 
(3-11) 

(3-12) 



In this way, the equation for f3 k is significantly simplified Mid reduces to equations 
similar to those we have solved in the initial value problem EJ) and in post-Newtonian 
studies. 113 Hereafter, we refer to the gauge condition presented here as an "approx- 
imate minimum distortion" (AMD) gauge. 

The above described AMD gauge condition is expected to have the following 
properties: 

• In a spherical symmetric spacetime with conformal flat initial conditions, it 
coincides with Smarr- York's MD gauge condition (Eq. (|3-2|) ) as well as with 
the MD gauge condition presented in this paper (Eq. ( |3-6[) ), because in this 
case, jij = 5ij and dtfij = 0. 



8 



• Using a post-Newtonian approximation as a guide, we can see that the difference 
between our AMD gauge and the MD gauge condition defined by Eq. (|3-6| ) 
appears, in general, at the third post-Newtonian order, because we neglect only 
coupling terms between hij (or and /^.E^'ll^ Hence, this difference is in 
general small, and we may expect that the global change rate of jij is still 
suppressed and is always sufficiently small. 

• As in the MD gauge conditions, in the AMD gauge condition, hij approximately 
satisfies a transverse relation (Fi = 0) everywhere except for a region near a 
rapidly rotating object as long as Fi is sufficiently small at t = (see Appendix 
B). Thus, we expect that such gauge condition will prove useful for analyzing 
gravitational waves in the wave zone, in which hij approximately satisfies a 
transverse-traceless (TT) condition, hijS 1 ^ = hij^& k = 0. 

In conclusion, the AMD gauge condition presented here possesses many of the merits 
of the MD gauge conditions, and it also results in equations which are much more 
easily solved numerically. 

§4. Analysis of gravitational waves 

In this paper, we analyze gravitational waveforms as measured by an observer 
located near the outer boundaries of the computational domain. Along the z-axis, 
this can be done through the quantities 

h+ = r(% x - % y )/2, (4-1) 
h x = r% y . (4-2) 

Since we adopt our AMD gauge condition and consider initial conditions with Fi = 0, 
hij is approximately TT in the wave zone. As a result, h + and h x are expected to 
be appropriate measures of gravitational waves emitted. (See §5 for a discussion 
of this.) The amplitude of gravitational waves will be largest along the z-axis as a 
result of our construction of initial data sets (see §5.1), and therefore h + and h x 
should be regarded as the maximum amplitude of gravitational waves. 

Gravitational waves are also measured through the gauge invariant Moncrief 
variables in flat spacetime, which are defined as follows:^ First, we transform from 
the Cartesian coordinates to the spherical polar coordinates, and then we split jij 
as 5ij + Y,im Clr> where & denotes 



Am 
Sij 



( H2Yi m h\i m Y[ m Q hu rn Yi mj(f 
r 2 (K lm Y lm + Gi m Wi m ) r 2 G im X lm 
\ * r 2 sin 2 6(Ki m Y lm - Gi m Wi m ) 



+ 



( -CimdyYim/ sm.6 C lm dgYi 

m sm a 

* r 2 D lm X lm /sme -r 2 D lm W lm sm8 | , (4-3) 

* * -r 2 AmA"/ m sin6* 



and * the relation of symmetry. f^Zm; hu m , K[ m , Gi m , Ci m and D\ m are functions 
of r and t, and are calculated by performing integrations over a two sphere of given 
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radius (see a previous paper llHP for detail). Yi m is the spherical harmonic function, 
and Wim and Xi m are defined as 



W lm 



(d g ) 2 - cot 6d e - — 27t(<9 ¥ 
sin 



Y, 



2& 



lm • 



dg — cot 9 

The Moncrief variables of even and odd parities are then defined as 



1 2(1-2)1 



(Z + 2)! 



4A^ m + Z(Z + l)fe 



l/m 



/ 2(Z + 2)! / q 
(Z-2)! I r 



where 



>4£m- 



K lm + Z(Z + l)G Jm + 2rd T Gi m - 2 
Hoi m 1 d 



llm 



2lm~~ 



2dr 



r{K lm + l(l + l)G lm } 



(4-4) 
(4-5) 

(4-6) 
(4-7) 

(4-8) 
(4-9) 



Note that using the above variables, the energy luminosity of gravitational waves 
can be calculated as 

dE _ r 2 
~dt ~ 32?r 

In this paper, we consider only the even parity Moncrief variables of I = \m\ 
which are expected to have the largest amplitude. 



{d t Rf m f + (d t R 



lm) 



(4-10) 



§5. Numerical study 

We have carried out numerical simulations of the merger between two identical 
relativistic clusters of collisionless particles for several sets of initial conditions. In 
the following, we present the results obtained and demonstrate that the AMD gauge 
condition discussed in §3 is robust and allows us to perform simulations for a time 
much longer than the dynamical time scale and to extract gravitational waveforms. 

5.1. Initial conditions 

For simplicity, we use conformally flat initial conditions. In this case, we only 
need to solve the Hamiltonian and momentum constraint equations, in addition to 
the equations for a and (3 k . The constraint equations are written as 



A 3 



-2itE^ - -AijA ij ilj- 7 
8 

87rJiV 6 , 

0) 



where = ifi Aij. As discussed in a previous paper, tiJ' we write Aij as 

2 

Ay = Wij + Wj fi - -SijW k:k , 



(5-1) 
(5-2) 

(5-3) 



10 



with 

W i = lB l -\( X ,i + B Kl x k \ (5-4) 



Then, the equations for Bi and x are derived from Eq. (5-2) as 

ABi = 8irJiTp 6 = 8irp*Ui, (5-5) 
A x = SirJiX^ 6 = -Mp*Uix\ (5-6) 



As a result of Eqs. ( p-g ) and (54), Eqs. (5-1), (5-5) and (5-6) are regarded as 
constraint equations to be solved. The strategy for setting up initial conditions is as 
follows: (1) we choose p*(x,y,z) and U{(x,y, z); (2) we solve Eqs. (5-5) and (5-6), 
from which Aij is subsequently obtained; (3) we solve Eq. (5-1). 

The rest mass density profile p*(x, y, z) of each cluster is simply chosen to be a 
parabolic function of the type 

/°o( 1 --r) r ± < r o, 



P*(r ± ) = { "V rl) '^' u ' (5-7) 
r± > r , 

where po an d r o are constants, and r± denote coordinate distances from the center 
of each cluster which is located at (±r c , 0, 0) with r c > r$. 

As one model, we choose the velocity field Ui for r± < tq as 

Ui (x, y, z) = uf b = (-O.Ol^F, ^V, o) , (5-8) 

where 



M one 

V = ^-£-, (5-9) 

and M° ne = 87rp r$/l5 is the total rest mass of one cluster. When we set the velocity 
field as Eq. ( |5-8[ ), we always refer to the model of the initial condition as Kepler 
model. In this model, we choose ujq = 1 (near Kepler angular velocity case) and 0.5. 
We also choose a nearly corotating velocity field as 

uAx 1 ) = uf h = f-0.01^y - —V, -V, o) . (5-10) 

V x r c r c J 

In this case, we set lvq = 1. We refer to the model of the initial condition as corotating 
model. 

Once p* and U{ are given, we can solve the constraint equations, as well as 
calculate the initial z-component of the total angular momentum of the system 

J = J d 3 xp*(xu y -yu x ) =j(oo). (5-11) 

Since we now consider the case in which each cluster is in near equilibrium, we 
must determine an appropriate velocity for each particle in addition to the orbital 
velocity. The method to determine the velocity (itj) for the a-th particle is as 
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follows. First, we numerically solve equations for an equilibrium state of a single 
spherical cluster assuming that is given by the parabolic function as Ecu (|5-7|) . 
The equilibrium state is easily obtained, as we have done in a previous paper. Once 
we obtain an equilibrium state of a spherical cluster, the square of the velocity field, 
u° ne u° ne = £ 2 , is obtained as a function of the coordinate distance from the center. 
Then, we assign the velocity of the a-th particle at (x a , y a , z a ) = (x p ± r c , y p , z p ) as 

(».,)„ H/>)( — cose, - , Vv sine) +u°J b , 

x l + vl ) 

( „„ )„ f ! j ( cos C + -== C) + 

, r p Jx p + y p Jx p + y p ) 




' y 2 

(u z ) a = -i(r P )^ -cosC, (5-12) 

r p 

where r 2 = x 2 p + yp 1 + z p , and random numbers of uniform deviation between and 
2tt are assigned for £. Namely, the velocity of each particle is written as a linear 
combination of the velocity obtained in constructing a single equilibrium spherical 
cluster and of the orbital velocity. In order to guarantee that the velocity field is 
given by Ui = u° rb at any point, we have distributed particles in an appropriate 
manner. Hence, in the limit r c — > oo, each cluster is in a true equilibrium. 

In addition to the constraint equations, we also need to solve equations for a and 
f3 k , where for the latter, we use the prescription discussed in the previous section. 
Since we require a to initially satisfy the maximal slice condition, K = = dtK, we 
solve the equation 

A{oajj) = 1va^(E + 2S k k ) + -a^ 7 A^A^ , (5-13) 

8 

after we have a solution for A%j and ip. 

5.2. Numerical results 

All simulations were carried out fixing parameters as ro = 1 and r c = l.lro, 
but using three different values for the total rest mass M* = 2Af° ne , ro/6, ro/4 and 
ro/3 (i.e., changing po). We note that in all cases, each spherical star is stable in 
isolation. Only when M* is greater than ~ ro/2.5, a spherical star becomes unstable 



in isolation for the density profile given by Eq. (5-7), because some orbits near the 
surface of the cluster are unstable. 

We typically use a uniform grid with (153, 77, 77) points and covering a co- 
ordinate domain with extents — 7.6ro < x < 7.6ro and < y,z < 7.6ro (i.e., 
5x = 5y = 5z = O.lro). In this case, radius of each cluster is covered by 10 grid 
points initially. We typically take the particle number N to be 10 5 Hereafter, all 

*' For some cases, we performed simulations changing N and the grid spacing, and found that 
the numerical results changed somewhat in particular after the collision of two clusters. This is due 
to the different magnitude of the fluctuations caused by the discreteness of the particle locations. 
However, the global features discussed here did not change considerably. 
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quantities are given in units vq = 1 (and G = 1 = c). 
5.2.1. Products of the merger 

In Figs. 1-3, we show snapshots of the particle positions at selected times for 
Kepler models with M^/tq = 1/6, 1/4 and 1/3, respectively. (We list several charac- 
teristic quantities of the initial conditions in Table I.) In the first two cases, the final 
configuration consists of a new nearly axisymmetric, rotating cluster, while in the 
third case a black hole appears to have formed. In the last case, the simulation could 
not be continued until we determined the apparent horizon because the accuracy be- 
came insufficient. However, we presume that a black hole was formed, because the 
merged object was very compact and a and 4> became, respectively, small (~ 0.1) 
and large (~ 0.7) when the simulation was terminated. The reason that we had 
to terminate the simulation before we determined the apparent horizon is probably 
related to our AMD gauge, which seems inappropriate when applied to a problem 
in which black hole formation proceeds and the forming region is not well resolved. 
We discuss this point further below. 

For the cases in which M*/r$ = 1/6 and 1/4, we stopped our simulations at 
t ~ 85ro when the new cluster seems to settle down to a nearly equilibrium state. For 
M*/ro = 1/6, the core within r < 2tq in the final state comprises about 82% of the 
total rest mass (i.e., m*(2ro)/M* ~ 0.82) and about 48% of the initial total angular 
momentum (i.e., j(2ro)/J ~ 0.48). For M^/tq = 1/4, we have m*(2ro)/M* ~ 0.83 
and j(2ro)/J ~ 0.61, respectively. It appears, therefore, that the final configuration 
consists of dense, rapidly rotating clusters surrounded by a halo. The reason that 
j(2ro)/J for M*/ro = 1/6 is smaller than that for M*/ro = 1/4 is probably related 
to the fact that in the first case, self-gravity is weaker, and particles of large angular 
momentum can easily diffuse outward. 

As shown in Figs. 1-3, the merging process toward the final state considerably 
depends on the initial compaction parameters. For example, for M^/ro = 1/6, 
1/4, and 1/3, the merging starts after about one, three quarter and half an orbital 
periods, respectively. Also, for M^/tq = 1/6, the merged object initially has a 
triaxial ellipsoidal shape which then gradually settles down to a spheroidal object 
forming spiral arms and transporting the angular momentum outward. On the other 
hand, for M*/ro = 1/4 and 1/3, the merging proceeds more quickly and violently, 
forming a very high density object or a black hole around the mass center. 

Since we did not prepare a quasi-equilibrium configuration as initial conditions 
for the mergers, it is difficult to present a rigorous interpretation of the results. How- 
ever, we believe that the different types of behavior at the initial stages of the merger 
can be related to one of the three following factors, (a) In all cases, each cluster in a 
binary is tidally elongated along the axis which connects the two cluster centers just 
after the simulation starts. As a result of the elongation, each cluster in the binary 
acquires a quadrupole moment, which makes the attractive force between the two 
clusters stronger and can trigger the merging. Since the elongation may be larger for 
binaries of larger compaction due to a general relativistic effect, the attractive force 
becomes stronger for clusters of larger initial compaction, (b) Gravitational radia- 
tion reaction, which accelerates the merging, should be stronger for more compact 
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binaries. This can be clearly deduced from the ratio ofjthe coalescence time to the 
orbital period, which can be approximately written asc3 



R T is ~ 1.4 for M^/ro = 1/3, while it is ~ 8 for M^/tq = 1/6. Thus, gravitational 
radiation can affect the evolution of the compact binary more efficiently, (c) For 
M*/ro = 1/3 and 1/4, the ratio of the initial orbital radius to the gravitational mass 
M is ~ 7 and 9, respectively, so that the orbital radius may be as small as the radius 
of the inner most stable orbit of the system. This implies that, for M*/r$ = 1/3 
and 1/4, a rapid transition from a stable near-circular orbit to a plunging orbit 
probably takes place when the orbital radius slightly decreases. On the other hand, 
for M^/ro = 1/6, this is never the case. 

Note that all the above factors are general relativistic. Therefore, although we 
cannot make a firm statement on the origin of the different types of behavior for 
different initial compaction, it is almost certain that a general relativistic effect, 
which is never observed in Newtonian simulations, is crucial in these simulations. 

The evolution in the late phase of the merger is also different for the three models. 
In particular, for M*/r$ = 1/4, we find an interesting feature which appears to be a 
peculiarity of the merger between collisionless clusters. In this case, the two clusters 
merge to form a very high density cluster in the early phase (t ~ 30ro). Then, a 
large fraction of particles expands outward (t ~ 30ro — 50ro), and the central density 
decreases. Because the system is strongly bound, the particles eventually converge 
in the central region (t > 50ro) to produce a final compact cluster. A possible 
explanation for the outward motion for t ~ 30 — 50ro can be found in the peculiar 
nature of a collisionless system, which does not possess any efficient mechanism of 
conversion of the kinetic energy into internal energy. In contrast, in the merger 
between fluid stars, a shock is formed, and the kinetic energy may be converted to 
the thermal energy. Consequently, most of the mass (or energy) in the system may 
be trapped near the mass center. 

For M*/ro = 1/4, we have also performed a simulation using the zero shift gauge 
condition {(3 k = 0). In this case, the large coordinate distortion leads to increasing 
large values of %j and produces very inaccurate results in a time ~ 1/3 orbital 
period. We show h xx and h yy in the equatorial plane for (3 k = and M*/ro = 1/4 
at selected times in Figs. 4(a) and (b), respectively. It is found that the maximum 
(minimum) value monotonically increases (decreases), and near t ~ 13ro it becomes 
~ 9 (—0.5) for h xx and ~ 2.4 (—0.5) for h yy . In contrast, in the simulations with our 
AMD gauge, the maximum (minimum) value of h xx and h yy does not monotonically 
increase (decrease) (see Figs. 4 (c) and (d)). Since we start the simulation with a 
conformal flat initial condition, it initially increases or decreases monotonically, but 
the maximum (minimum) stops increasing (decreasing) in about one orbital period 
when they have reached ~ ±0.05. This behavior is a clear indication that the AMD 
gauge condition provides the required suppression of the coordinate distortion. 

Note that the magnitudes of htj are roughly consistent with h{j being a second 




(5-14) 
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post-Newtonian quantity, i.e., of 0[(M*/ro) 2 ]. This confirms that our AMD gauge 
condition is actually a valid approximation of the MD gauge conditions, because the 
coupling terms between hij and /3 fc are sufficiently small. Since hij does not mono- 
tonically increase nor decrease and the absolute value of each component remains 
small, simulations can be stably and accurately performed for a time much longer 
than the dynamical time scale using our AMD gauge condition. 

In Figs. 5 and 6, we show snapshots of the particle positions at selected times 
for corotating models of M*/ro = 1/4 and 1/3. These figures should be compared 
with Figs. 2 and 3, respectively. As in the Kepler case, the final product of the less 
compact binary is a rotating cluster, while that of the more compact binary appears 
to be a black hole. In the case M*/ro = 1/3, again the simulation could not be 
continued until the apparent horizon was determined. 

By comparing Figs. 2 and 5, we find that the merging process in corotating 
models is different from that occurring in Kepler models, and this is particularly 
evident for the case M^/tq = 1/4. In general, in the cases of corotating initial 
conditions, the outer part of the binary has more angular momentum than that 
for the Kepler model, so that the centrifugal force is stronger. Hence, the merging 
proceeds mainly in the inner region of the system, and particles in the outer region 
are not involved. Since the two clusters do not collide very quickly and the merging 
proceeds gradually, the resulting object does not expand outward significantly after 
the merger when compared with the Kepler case. The product of the merger settles 
down to a spheroidal object more gradually, forming spiral arms in the outer region. 
In this case, m*(2ro)/M* ~ 0.79 and j(2ro)/J — 0.48 in the final stage at which we 
terminated computation, i.e., at t ~ 80ro- The fact that the final value of j(2ro)/J 
is smaller than that for the Kepler model with the same initial compaction results 
from the initial corotating velocity field, in which the outer parts have a large angular 
momentum and are therefore able to diffuse outward more easily. 

In order to clarify the reason that both our AMD gauge and the MD gauge con- 
ditions are not well suited to study the formation of a black hole, we also performed 
simulations for a Kepler model of loo = 0.5 with M*/ro = 1/3. In this case, the an- 
gular momentum is not large enough to counteract the gravitational force between 
the two clusters, and they quickly merge into a black hole. In this simulation, we 
set N = 50000, and the grid spacing was 5x = ro/15. 

As in the Kepler and corotating cases with loq = 1 and M*/ro = 1/3, the 
simulation using the AMD gauge condition terminated before an apparent horizon 
was formed. The reason for this seems to be the following: Using a MD gauge 
condition or our AMD gauge condition with K = 0, Eq. (3-7) can be roughly 
approximated as 

Ap\ i ~ 12vrJ M . (5-15) 
During a gravitational collapse we have 

J r < 0, and J M ~ ^d r (J r r 2 ) < 0, (5-16) 

where J r denotes the radial component of Jj and we assume that the matter collapses 
mainly in the radial direction. Thus, in the central region where the black hole 
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formation takes place, the divergence of /3 l is positive, i.e., 

(3\i > and (f > 0. (5-17) 

Clearly, the condition (3 r > implies that coordinates spread radially outward; 
i.e., the physical size between two neighboring grids increases and the numerical 
resolution deteriorates. At r = 0, <fi increases rapidly when the collapse begins 
(see the solid line in Fig. 7); for t ~ 14ro, 4> ~ 0.8 while it was ~ 0.13 at t = 0. As 
a result of this variation in (j), the physical separation between two neighboring grid 
points around r = increases by a factor of e ~ 4. If we require the same accuracy 
as that at t = 0, we have to reduce the grid spacing 5x by a factor of 4, which is not 
feasible in a 3D numerical simulation. This problem might be resolved with the help 
of an adaptive mesh refinement (AMR) technique. 0) However, (f> increases more and 
more as the collapse proceeds, so that even an AMR treatment might be insufficient 
without changing the gauge condition once a gravitational collapse begins. In this 
paper, thus, we propose a modified spatial gauge condition derived from our AMD 
condition as 

k 

P h = P k AMD-y(3 r AMDf(.r,t), (5-18) 

where P^md 1S the shift vector determined from our AMD gauge condition, (5 r = 
x fc /3^ MD /r, and f(r,t) is a function of r which satisfies f(r = 0) = O(l) and 
f(r = oo) = 0. In this case, (3 r around a black hole forming region can be set 
to a small value. Although the distortion in 7^ due to the radial motion increases, 
the coordinate twisting due to the angular motion may still be suppressed. Fur- 
thermore, the TT property is preserved in the wave zone because f(r,t) — > for 
r — > 00. To validate this idea, we perform a simulation using the following f(r,t) as 
an example, 

f( r ,t) = f (t) . ] , r XR , (5-19) 

J v ' ; ,M,W 1+ (r/3M*) 6 ' v ' 



where we choose /o to be 



fo(t) = \l f T a{ J = $ - °t (5-20) 
J w 10 for a(r = 0) > 0.5. v ' 



In Fig. 7, we show <fi at r = for /o of Eq. (5-20) (the dashed line) as well 
as for /o = (the solid line). The behavior of <p clearly shows that the coordinate 
spreading is suppressed after we turn on a non-zero /o- In Fig. 8, we also show the 
root mean square of the coordinate position for particles defined as 



\ 



^E«- (5-21) 

i,a 



Although the physical results of these simulations should be identical, the time evo- 
lution of r rms in each simulation should be different when different gauge conditions 
are adopted. Since the coordinate position of each particle becomes smaller when 
the radial coordinate spreads outward, we expect r rms to be smaller for a smaller 
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/o- Figure 8 clearly reflects this property and that the non-zero /o suppresses the 
coordinate spreading effect. 

In Figs. 9 and 10, we also show snapshots of particle positions at selected times 
for /o = and the non-zero /o, respectively. We find that in the fo = case and at 
late time, particles are concentrated in very narrow regions around the origin and, 
as a consequence, the resolution becomes insufficient for t ~ IAtq. On the other 
hand, in the non-zero fo case, the coordinate spreading was not appreciable, at 
least up to the formation of an apparent horizon, and the simulation continued after 
the formation of the apparent horizon and for a long enough time that the horizon 
stretching became dominant and spoiled the numerical accuracy. In conclusion, we 
can mention that in addition to our AMD gauge condition, it may be helpful to 
incorporate an appropriate non-zero function f(r, t) in the final phase of the merger 
just before the formation of a black hole. 

It is important to stress that the present strategy to compensate for the limita- 
tion of MD-type gauge conditions is well suited to the problem under investigation 
and is not the only possible solution. In a different scenario, different solutions might 
turn out to be more effective. Consider, for example, the case in which a rapidly 
rotating object collapses and a disk is formed prior to the formation of a black hole. 
In this case, it would probably be more appropriate to suppress /3 Z first, and then, 
f3 R [= (xP x +y(3 y ) j \J x 1 + y 2 ] only when the black hole formation starts. It is evident 
that more detailed studies on the choice of the gauge condition suitable to study the 
formation of a black hole are necessary. 

5.2.2. Gravitational waveforms 

In Figs. 11 and 12, we show h + /(M 2 /r ) and h x /(M 2 /r ) at z Q bs = 7.5r (solid 
lines) and 6ro (dashed lines) for the Kepler and corotating models with M* = ro/4, 
respectively. We plot these lines as a function of (t — 2 o bs)/ r 0- Hence, if h+ and 
h x behave as solutions of the wave equations near the outer boundaries, the two 
lines should approximately agree. Note that gravitational waves begin to reach the 
observer at t — z Q bs — 0, so that the waveforms shown for t — z b s < are meaningless 
as the coalescence waveforms. 

It should be mentioned that the expected wavelength of gravitational waves in 
the initial phase of the merger is approximately T/2 ~ 20ro, where T is an approx- 
imate orbital period (see Table I for definition of T). In principle, the asymptotic 
waveforms should be extracted in the wave zone, and therefore for z Q b s > 20ro. Nev- 
ertheless, the waveforms shown here appear to constitute a fair description of the 
asymptotic waveforms due to the following reasons: (i) The solid and dashed lines 
agree well in both cases throughout the computations, so that h + and h x actually 
propagate at the speed of light, (ii) The wavelength in the initial phase of the 
merger approximately agrees with that expected from the initial orbital period, (m) 
The amplitude of h + and h x in the initial phase roughly agrees with an analytical 
estimate using the quadrimole formula for two point masses in a circular orbit, i.e., 
M*(47rr c /T) 2 = M 2 /2r c .i$ (iv) The waveforms agree fairly well with those expected 
from Figs. 2 and 5; i.e, in their initial merging phases, the amplitudes of gravita- 
tional waves are the largest because the systems are highly non-axisymmetric, but 
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once they settle down to form an elliptical merged object, the amplitude gradually 
decreases. 

The gravitational waveforms in Figs. 11 and 12 also reflect the difference in the 
merger process for the Kepler and the corotating models. In the Kepler case, the 
merging quickly proceeds to form a new, high density rotating object resulting in a 
rapid and large increase in the amplitude of gravitational waves. In the corotating 
model, on the other hand, the merging proceeds less rapidly, and the resulting am- 
plitude of gravitational waves is smaller. Also, in the Kepler case, a high density and 
rapidly rotating ellipsoidal object is formed after the merger, and it is responsible 
for the emission of high frequency gravitational waves. The results presented here 
suggest that waveforms of gravitational waves from the merger of stellar clusters 
may be sensitive to the internal velocity fields that two clusters possess before the 
merger, and could be used to extract astrophysical information. 

In both the Kepler and the corotating cases, the maximum values of h+ and 
h x are roughly ~ 0.08M*(4M*/ro) and the wavelength is ~ 10ro(ro/4M*) 1 / 2 , as 
expected from the quadrupole formula. If similar mergers between highly relativistic 
clusters actually occurred in the early universe, the amplitude of gravitational waves 
could be 

-- 18 (^)fe)(^)- 

with the frequency ~ lO~ 2 (lO 6 M /M*)(4M*/r o ) 3 / 2 Hz, where M denotes the solar 
mass. If mergers between highly relativistic star clusters of mass ~ 1O 6 M and 
radius ~ 10M* occurred in the early universe, it should be possible to detect emitted 
gravitational waves by the planned gravitational wave detectors in space. U> 

In Figs. 13 and 14, we show ri?2,2±/ r o as a function of (t — r b s )/?"o f° r the 
Kepler and the corotating models of M* = ro/4, where we define 

R2,2 + ^ R ^ + J>>-\ (5-23) 

R2,2- = R ^~^-\ (5-24) 

The solid and dashed lines denote the waveforms extracted at r Q b s = 6.5r and 7.3r , 
respectively. 

We find that for ri?2,2±) the solid and dashed lines do not coincide well compared 
with those for h + and h x . In particular, for t — r Q b s < 10ro, in which case the 
variation time scale of the system is relatively long, the coincidence is not good. The 
disagreement is probably due to the fact that the main contribution to i?2,2± comes 
not only from the 0(r _1 ) part of 7™ but also from the 0(r -3 ) part of exp(</>), which 
can be approximated for large r as 

exp(</>) = l + y r + ^fj^ + 0(r- 5 ), (5-25) 
where denotes the trace free part of a quadrupole moment. An order of magni- 
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tude estimate for / = \m\ = 2 modes of hij and exp(0) is given by 

hij ~ 7~dt^ 1, and exp( ^ ~ 73- (5 ' 26) 

As a result, the ratio of the magnitude of hij to that of exp(</>) is ~ (wr) 2 , where 
u) denotes a characteristic angular frequency of gravitational waves. In order to 
extract gravitational waves (i.e., the contribution from hij) clearly, tor should be 
much larger than unity. However, for t — r ^ s ~ 0, we extract R% t 2± for r ~ 
and the magnitudes of the two contributions are roughly equal. For this reason, the 
solid and dashed lines do not coincide well, and R,2,2± for t — r Q b s < 10r should not 
be regarded as accurate quantities for measuring gravitational waves. 

However, when u is sufficiently large, the coincidence between the solid and 
dashed lines is improved, and for t — r b s ~ 10 — 40ro, the coincidence is fairly 
good. Assuming the quantities R.2,2± describe gravitational waveforms correctly, we 
can calculate the total energy luminosity AE due to I = \m\ = 2 components by 
integrating dE/dt from t— r Q b s = 10r to the end of the simulations using Eq. fl4T0]) . 
We obtain for the Kepler and corotating models, respectively, 

AE ~ 4.5 x 1(T 3 and 1.5 x 10 _3 M*. (5-27) 

For a binary composed of equal masses in a circular orbit of orbital radius 2r c , the 
energy radiated by gravitational waves in one orbital period is roughly given by the 
quadrupole formula as 

which confirms the validity of the results obtained in Eq. ( |5-27| ). The reason that 
AE for the Kepler model is larger than that for the corotating model is that the 
merging proceeds more quickly and that a rapidly rotating ellipsoidal core of a higher 
density is formed in the initial stage of the merger. Irrespective of the used model, 
however, these results apparently indicate that a merger between two relativistic 
clusters is an efficient source of gravitational waves. 



§6. Summary 



We have performed fully general relativistic simulations of the merger between 
relativistic clusters of collisionless particles adopting a new spatial gauge, the AMD 
(approximate minimum distortion) gauge condition. We have reached the following 
conclusions: 

• By using the AMD gauge condition, it is possible to perform numerical simu- 
lations of coalescing binary clusters for a time much longer than the dynami- 
cal time scale of the system when the merger does not produce a black hole. 
Also, it is found that with this gauge condition, gravitational waveforms can be 
calculated fairly accurately. These results suggest that with the AMD gauge 
condition, the spatial coordinate distortion is suppressed to a level adequate to 
carry out stable and accurate simulations. 
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• Our AMD gauge condition, as well as a class of spatial gauge conditions similar 
to the MD gauge conditions, still has an undesirable, coordinate spreading 
property (i.e., (3 r > 0) when applying it to the study of black hole formation. 
As a result of this property, the spatial resolution around a black hole forming 
region rapidly tends to become insufficient. 

• However, if our AMD gauge condition is slightly modified during the process 
of black hole formation, it is still possible to carry out simulations up to the 
formation of an apparent horizon. Using this approach, we have shown that 
it is possible to perform numerical simulations over a sufficient long time scale 
from the early merging phase up to the formation of a black hole. 

• In a merger between two highly relativistic star clusters with M*/ro = 1/4 into 
a new cluster, the amplitude of gravitational waves has been estimated as ~ 
lCT 18 (A-iyiO 6 M ) at adistance of ~ 4000Mpc with a frequency ~ 10" 2 (10 6 M Q /M)Hz. 
Also, ~ 0.5% of the rest mass energy may be dissipated by gravitational waves 

in the final phase of the coalescence. 

The spatial gauge conditions presented and investigated in this paper are useful 
for a wide variety of problems, and we expect that they will be widely exploited 
for performing simulations of coalescing binary neutron stars, which represent the 
most promising sources of gravitational waves for kilometer-size laser interferometers. 
When a general relativistic hydrodynamic code is installed to replace numerical code 
for collisionless particles, we will be able to perform the numerical simulations for 
such binary systems and follow them up to the formation of black holes or new 
rotating neutron stars. This work is now in progress, and in future papers we will 
present the numerical results. 123 
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Appendix A 

Numerical treatment for the transport terms 



For treating the transport terms in the evolution equations ofgeometric vari- 
ables, we use a method similar to that adopted by Stark and Piran.E3) In the follow- 
ing, we demonstrate the treatment only in the x direction, but the same operations 
are also carried out in the y and z directions. We also assume a uniform grid along 
the x direction and denote the grid spacing and time step as 5x and 5t. 

Evolution equations for the geometric variables may be written in the form 

(d t -p x d x )Q = S, (A-l) 
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where Q denotes one of the geometric variables and S denotes the source term. The 
solution for the finite difference equation is written 



Q 



n+l 

i 



HQ?+i - QU) - W\(QU - 2Q? + Q?) 



in the first order upwind scheme, and 



Q 



n+l 



Q? 



1 



HQU - QU) - »\QU - 2Q™ + Qti) 



+ S? +1/2 5t, (A-2) 



+ S? +1/2 6t, (A-3) 



in the second order scheme. Here, v = —{(5 x ) n t +l ^ 2 5t/5x, and Qf denotes Q at the 
n-th time step and i-th grid point. In this paper, we compute Q™ +1 as 



Q 



n+l 



v{QU - Q?_x)-(1 - *)M(G?+i - 2Q? + Q?_ 



-s l v 2 (Q? +1 -2Q? + Qt 1 ) 



+ +1/2 <ft.(A-4) 



Here, Sj is a limiter function at the i-th grid point, chosen as 

25Qj6Qi-i + e 
5Qf + 5Q|_i + e 

for 8Qi6Qi-i < 0, 

where 5Qj = Q™ +1 — Q™ and e is an appropriately small constant. 

Appendix B 

Behavior of F{ in our AMD gauge 



(A-5) 



Substituting Eq. hito Eq. (2-20), we obtain 

(d t - (3 k d k )F l= 2a(f kj A ik)j + f kj tj Ai k - Ji^ iM ) + 2f k a k A. 



+6 jl p" il (h ij , k + h lkJ ) + h u p\ jk V k + Fi(3 l ti + hflp^P 1 



= sr. 



(B-l) 



Here is a nonlinear function and can be set to zero in the linear approximation. 
Hence, in this approximation, Fi is always zero if it is zero initially. Even when we 
consider a linear perturbation in the Schwarzschild spacetime or the spacetime of a 
spherical star, Sf is non-linear and F{ is always zero if initially zero. However, in a 
spacetime of a rotating black hole or star, hij (/*-'), (i k and Ajj appear at zeroth order, 
so Fi also appears from zeroth order. Nevertheless, in the case that the rotation of 
the object is not very rapid, hy (3 k and Ajj are small, and hence F{ can be 

considered to be small. Even in the case that the object is rapidly rotating, hij 
[3 k and Ajj quickly decay as 0(r~ 2 ), 0(r~ 2 ) and 0(r -3 ), respectively, so that Sf 



21 



and Fi can still be regarded as small quantities everywhere except near the rapidly 
rotating object. Furthermore, in the post-Newtonian approximation, Sf contains at 
most third post-Newtonian quantities, because hij (f 13 ), f3 l , and Aij canbe regarded 
as the second, first, and first post-Newtonian quantities, respectively. Thus, 
except in the highly relativistic region in which hij, f3 l and become O(l) or a 
relativistic region near the rapidly rotating relativistic objects, Sf is guaranteed to 
be small, and hence, Fi is expected to be small as long as it is small enough on 
an initial time slice (i.e., as long as we do not consider strange initial conditions in 
which Fi is large initially). 
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Table I. The list of initial conditions and final states for simulations performed 
in §5. M* = 2M™°, M, J and T denote the total rest mass, gravitational mass, 
total angular momentum and approximate initial orbital period (27ry / 8rJ/M*), re- 
spectively. All the quantities are given in units r$ = I (and G = 1 = c). In the cases 
marked with f , we could not determine the apparent horizon formation, but black 
holes seem to be formed. 





M/M, 


J/M 2 


T 


velocity field 


Final state 


Figures 


1/6 


0.966 


0.974 


50.2 


Kepler (loo = 1) 


rotating cluster 


1 


1/4 


0.952 


0.818 


41.0 


Kepler (u>o — 1) 


rotating cluster 


2 ,4, 11, 13 


1/3 


0.941 


0.726 


35.5 


Kepler (u>o — 1) 


black hole T 


3 


1/4 


0.955 


1.01 


41.0 


Corotation 


rotating cluster 


5, 12, 14 


1/3 


0.944 


0.894 


35.5 


Corotation 


black liole 1 ^ 


6 


1/3 


0.931 


0.371 




Kepler(cJo = 0.5) 


black hole 
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